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Abstract 

Continuum strain energy density functions are developed for soft biological tissues that 
possess slender, fibrillar components. The treatment is based on the model of an elastica, which 
is our fine scale model, and is homogenized in a simple fashion to obtain a continuum strain 
energy density function. Notably, we avoid solving the exact, fourth-order, nonlinear, partial 
differential equation for deformation of the elastica by resorting to other assumptions, kinematic 
and energetic, on the response of individual, elastica-likc fibrils. The formulation, discussion of 
responses of different models and comparison with experiment are presented. 



1 Background 



Currently- used strain energy density functions for soft biological tissue have two main origins. 
Some have been adopted from the rubber elasticity and polymer elasticity literature, and oth- 
ers have functional forms tha t have been chosen to reproduce the characteristic locking behavior 
observed in experiments (see Fung . 19931 . for a detailed treatment). Among the rubber/polymer 
elasticit y models are micromechan ically-derived ones, which mostly incorporate entropic elastic- 
ity (see Landau and Lifshitzl . 1951 . for a discussion). Entropic-elasticity models are suitable for 
materials in which the uncoiling of long chain molecules under axial force causes a decrease in 
configurational entropy as few e r con figurations become available to the molecule vibrating under 
its thermal energy (see lOgdenl . 119971 . for detailed treatments of these models). However, it is not 
clear that the application of entropic elasticity is appropriate for many soft biological tissues such 
as tendons, ligam ents and m u scles. As an example, consider the case of tendons, which have a high 
collagen content. Sun et all (120021 1 demonstrated by laser trap experiments that the elasticity of 
the collagen molecule, which is a triple helix with a diameter of 1.5 nm and a contour length (fully 
uncoiled length) of appro ximately 300 nm, is well-represented by the Worm-like Chain Model of 
Kratkv and Porodl (|l949l b However, collagen is not restricted to the form of long chain molecules 
in the tendon. It forms fibrils of around 300 nm diameter, and lengths of the order of 100s of /mi. 
These are further ordered into fibers that can run the entire length of tendons (the order of cm). 
The entire hierarchical structure has extensive crosslinking, including a longitudinal staggering of 
the collagen molecules that leads to a chara cteristic band e d structure on the scale of a fibril , and a 
"crimp" with a wavelength of 10 — 50 fim. Screen et all » IProvenzano and Vande^ B. 
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Given this ordered, hierarchical structure with extensive crosslinking, one must question the use of 
entropic elasticity. Due to kinematic constraints imposed by the crosslinking it seems unlikely that 
the collagen molecules are able to sample many configurations via thermal fluctuations. Similar 
arguments can be made for ligaments and muscles. 

n (|l987T l. Strain- 



Support for this view may be inferred from the experiments of IWoo et al. 



controlled cyclic tension tests of canine medial collateral ligaments at temperatures between 2° C 
and 37° C showed that the area of hysteresis loops on the stress-strain diagram decreased as the 
temperature of the experiment increased. The use of strain control implies that the decrease in 
hysteresis was associated with a reduction in initial modulus of the stress-strain curve. In a standard 
viscoelastic solid, the initial modulus is a material property; in particular, it is independent of 
viscosity, and therefore not subject to mechanisms of relaxation that may be perceived as a decrease 
in modulus. This argument leads to the conclusion that, in these experiments, the initial modulus 
of the stress-strain response was decreasing with an increase in temperature. A decrease in elasticity 
with increase in temperature is a signature of elastic modulus that arises from variation of internal 
energy, not of entropic elasticity, which makes the initial modulus increase with temperature (see 



Treload . fll)75T ). 



One reason for the attractiveness of entropic elasticity models is that they reproduce the 
experimentally-observed response of soft biological tissue in uniaxial tension shown in Figure 
We draw attention to the characteristics: a prolonged initial regime with low modulus (the "toe" 
region), followed by a short nonlinear regime with rapidly- increasing modulus (the "heel" region) 
and a final high modulus region. We will refer to this as the "characteristic soft tissue response". 
However, there are other, non-entropic, models that also reproduce this behavior. It has been typ- 
ical in the biomechanics literature to use strain energy density fun ctions with m athematical forms 
that are designed solely to possess this characteristic response (see iFund . 119931 ) . This second class 
of models is, however, limited by the lack of microstructural bases for the corresponding strain 
energy functions. 

This paper is founded on the recognition that the characteristic locking behavior can be modelled 
by an internal energy-based, i.e. non-entropic, model of elasticity that accounts for the uncoiling 
of crimped fibrils with increasing tension. In these models, the characteristic soft tissue response 
is therefore determined by the elastica-like force versus tip displacement response of the individual 
fibrils. This consideration leads to a micromechanically-derived strain energy density function for 
soft tissue, which, as argued above, has the proper basis in internal energy, rather than entropy 
effects, which are suppressed due to crosslinking. 

The realization that characteristic soft tissue response is mode lled by the forc e -displ acement 
response of an elastica — or approximations of it — is hardly new. Diamant et al. ( 1972! ) used a 
planar model of rigid links joined by elastic hinges, which they related to the el astica, to model 
their observations of stress-stretch behavior of rat tail tendons. In [Dale et al.l (| 19721 ) four kinematic 
models of crimped fibers were considered: a planar sinusoidal waveform, a helical shape, a zig- 
zag waveform with hinged apices and a zig-zag with apices that undergo bending to maintain a 
constant angle whi le deforming. The change in profile of these waveforms was studied and compared 
with experiment. Beskos and Jenkins ( 19751 ) modelled mammalian tendon as an incompressible 
composite with a continuous distribution of inextensible fibers with a helical shape. The assumption 
of inextensibility dominates the response of this model leading to s tress locking in uniaxial tensio n 
at a finite stretch. The planar assumption was also adopted by Comninou and Yannasl ( 19761 ). 
who modelled single collagen fibers as sinusoidal beams. Using the theory of shear deformable 
beams with linear constitutive relations for axial stretching and bending, but allowing geometric 
nonlinearities, they obtained a nonlinear stress-strain response of single fibers and extended it to a 
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composite with uniaxial reinforcement by sinusoidal fibers. In Lanir ( 19781 ) a planar model of a beam 
on an elastic foundation was adopted for the mechanical interaction between collagen (modelled a s 
a beam) and elastin (the elastic foundation). Similar ideas were explored by Kastelic et al. (1980), 
in whose model crimped collagen fibers were modelled by links that have neg li gible stiffness until 
fully extended. The classical theory of elasticas was used by iBucklev et all (11 980) to treat the 
deformation of slender filaments, and the model was solved numerically. iBasu and Lardner (|l985T ) 
also studied the stress-stretch response of sinusoidal beams in elastic matrices, although they did not 
make the li nk to fibrous soft tissu e. A kinematic chain with finite axial stiffness and torsional springs 
was used bv lStouffer et al.l (119851') to represent the u ncoiling of crimped fibers, and compared against 
experiments. More recently, iHurschler et al.l (|1997l ) developed a strain energy density function for 
tendon and ligament with seven parameters including microstructural organization to describe 
the stress-stretch behavior. The authors also derived simplified versions of t heir model that were 
used t o fit experimentally-determined, nonlinear stress-stretch curves. Finally, iFreed and Doehring 
(j2005r i returned to the assumption of a helical structure for collagen fibrils, and using Castigliano's 
theorem, obtained the force-displacement relationship. 

In this communication, we present a very general and powerful procedure for developing the 
strain energy density function for soft tissue based on the elastica as a model for slender fibrillar 
structures. To fix ideas, we refer to collagen fibrils. In a notable departure from the body of 
work cited above, we first obtain the exact, nonlinear, fourth-order elliptic partial differential 
equation for the quasi-static deformation of the extensible elastica. The underlying kinematics are 
fully nonlinear and the elastica's strain energy is assumed to be given by quadratic functions of 
curvature and the Green-Lagrange strain tensor0 The difficulty of obtaining analytic solutions of 
even simpler versions of the governing partial differential equation has been noted by some of the 
authors cited above. Furthermore, numerical solutions, while possible, will prove both expensive 
and cumbersome, since our ultimate aim is a strain energy density function for composite soft tissue 
in which the elastica-like fibrils are the reinforcements at a microscopic scale. For this reason we 
have examined a few distinct assumptions that make it possible to obtain force-extension solutions. 
These assumptions are related to the kinematic and energetic behavior of the microscopic collagen 
fibrils. At the outset we wish to emphasize that these assumptions, in addition to being well- 
motivated in our view, are most important to this study because they deliver a tractable problem, 
with only a small number of parameters that are also physically-meaningful. This is a theme that 
we will return to at several points in the paper. With this approach we have also been able to 
identify a model that can be made to correspond well with experimental data. A rigorous validation 
of these assumptions must, however, await in situ studies of deforming fibrils. 

The organization of the remainder of the paper is as follows: In Section [2] we lay down the 
fundamental problem of the elastica. The cases of elasticas that are restricted to circular and 
sinusoidal arcs, and have further constraints of kinematics and energetics imposed upon them, are 
developed in Sections [3] and HJ respectively. The various models for the deforming elastica are 
compared against each other, and against experiment, in Section [5l The extension to macroscopic 
strain energy density functions, from which the tissue stress-stretch response can be obtained, and 
a basic discussion on convexity appear in Section El Closing remarks are made in Section 



lr The latter dependence is motivated by the St. Venant-Kirchhoff model, but there is no further significance to 
this choice. Others are equally admissible and do not imply substantive changes in the outcome. 
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2 The deforming elastica 



Consider the elastica, a curve, T C M 3 , parametrized by its arc length coordinate, S, in the reference 
configuration (Figure [T|). Points along the curve are identified by position vectors X(S) £ M 3 . The 
tangent at S is T(S) = dX/dS. Using a Cartesian basis of orthonormal vectors {ei, e2, 63} it is 
clear that dXj/dS is a direction cosine, say cos a/, where I = 1,2,3, and ai are the corresponding 

angles of inclination of T. Using the Euclidean norm of a vector, ||t>|| = \J^2i vf, it then follows 

that ||T|| = 1. The curvature of r is kq = \\d 2 X /dS 2 \\. In the deformed configuration, 7, points 
have position vectors x(S) = X(S) + u(S). The tangent vector to T is carried to t(S) = dx/dS, 
and by the chain rule it can be written as 

dx dX 




Figure 1: The curve, T G M 3 . 

Introducing the deformation gradient tensor, F := dx/dX , gives t = FT, from which it follows 
that ||t|| 2 = T ■ CT, where C = F T F is the right Cauchy-Green tensor. The stretch along S will 
be denoted by A := ds/d5 = ||t||, and satisfies A > 0. As an aside, note that if the arc length in 
the deformed configuration, say s, is used to parameterize x as x(s), then the tangent with respect 
to 7, which we will denote by t*, is given by = da;/ds, and satisfies = 1. The curvature 
of 7 is k = ||d 2 s/ds 2 ||, and, clearly, it can be written as k = ||d 2 a;/A 2 dS' 2 ||. 

Motivated by the decoupling of the bending and stretching energies of the classical Euler- 
Bernoulli beam, we assume that in the present nonlinear setting the strain energy of T can also 
be decomposed into bending and stretching components, W(k, A) = K(k) + U(X). The bending 
stiffness is denoted by B, the stretching modulus by E, cross-sectional area by A, and the reference 
contour length by L. This gives 

L L 

K(k) = j ±B(k - Ko ) 2 dS, U(X) = j \eA(\ 2 - ifdS. (1) 
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Observe that W(k, A), K(n) and U(X) denote functionals of the corresponding arguments. For 
simplicity, we will assume that the bending and axial stiffnesses are constant; B = const., EA = 
const. Also see Section [2.11 in this regard. 

The governing partial differential equation for the deformation of the elastica is obtained by 
imposing stationarity of the following free energy functional, in which the above definitions of 
curvature and stretch have been exploited: 



u 




d 2 



X 



ds 2 



Kq(S) 



+ -EA 
2 



dx 



dS 



q ■ u dS. 



(2) 



Note that the strain energy has been incorporated from (TjQ) and q(S) is the external force per unit 
contour length. We have also assumed that the displacement is specified at S = {0, L}: u(0) = 
and u{L) = g. Introducing the variations x e = X + u + ew, where e € M. and w(S) is an arbitrary 
vector field satisfying w(S) = at S = {0, L}, the stationarity condition is 



de 



Sf[u 6 



e=o de 



-B 



d 2 x, 



ds< 



ko(5)) +-EA 



dx F 



dS 



1 | -q u £ | dS I - U. 

£=0 



(3) 



Standard variational calculus, the arbitrariness of w and dw/dS, as well as the homogeneity of 
w at S = {0,L}, yield the following Euler-Lagrange equations: 



B- 



Xd(XdS) 
with the boundary conditions 



(1 - Kq/k) 



d 2 X \ 



Ad(Ad5) 



(4) 



u(0) = 0, u{L) =g, B ( 1 



d 2 a; 



k J Xd(XdS) 



at 5 = {0,L}. 



(5) 



Since A = ||cLc/d/S r || and k = ||d 2 a;/A 2 dS' 2 ||, Equation ^ possesses complexity beyond that 
apparent in its form above. In addition to the boundary conditions on u at {0, Lj, note that 
the generalized force satisfies homogeneous boundary conditions at S = {0,L} in ((5])3 o Alternate 
boundary conditions can be used in the above variational procedure. 



2.1 Simplifying assumptions on kinematics and energetic response of the elas- 
tica 

The ultimate aim of this work is the development of a macroscopic strain energy density function, 
where the micromechanics arises from the deformation of the elastica. The exact micromechanics 
of the deforming elastica is obtained by solving the partial differential equation subject to the 
boundary conditions in ([5]). Its solution, however, is nontrivial on account of the nonlinearity and 

2 This generalized force is conjugate to dw/AS in the Euler-Lagrange equations arising from J3J, and therefore 
admits the interpretation of a moment. 



5 



fourth-order form of the partial differential equation. It is desirable to avoid the complexity and 
expense of solving this equation repeatedly in fine scale computations that will be coarse-grained 
in some suitable fashion to obtain the strain energy density function and stress at each material 
point on the macroscopic scale. For this reason, we examine certain assumptions, kinematic and 
energetic, which simplify the micromechanics. 

The kinematic assumption in Section [3] is of an elastica which deforms through a family of cir- 
cular arcs. In Section 0] it is a family of sinusoidal waveforms. In both sections we also consider two 
further kinematic assumptions, inextensibility and planar incompressibility, and the assumption of 
a stationary strain energy. Planar waveforms h ave been used in ma ny of the studies cited in Section 



D and reported in the experimental studies of I Screen et al.l ( 20041 ) and lProvenzano and Vanderbv 



200 6|) . Helical forms ha ve also been reported in a few cases (see iBeskos and Jenkins! . Il975l ; 



Freed and DoehrineJ . 20051 . although even these authors cite at least an equal number of papers 



that reported planar waveforms). Since there exists some dispute in the literature in this regard, 
we have chosen the simplicity of planar waveforms. We also note the ease of parametrization of the 
circular and sinusoidal forms. 

The persistence of the initial family of waveforms, and the assumptions of inextensibility or 
planar incompressibility are, perhaps, the strongest assumptions in this paper. The assumption on 
persistence of the family of waveforms maintains the ease of parametrization gained by assuming 
circular arc and sinusoidal initial waveforms. For each of the assumed waveform families (circular 
arc and sinusoidal) and additional kinematic assumptions (inextensibility and planar incompress- 
ibility) we treat the effective elastic parameters, B and E, as those measured in an experiment in 
which the elastica's deformation remains within the waveform family, with the additional kinematic 
assumptions holding. We note that the effective elastic parameters so measured will be different 
from those obtained in unconstrained experiments. For simplicity we assume B and E to be con- 
stant for each combination of waveform family and kinematic assumption. Also note, however, 
that these parameters will differ in each one of these cases. In each case, we will not solve the 
unconstrained problem, but will restrict ourselves to the specific waveform family, and additional 
kinematic assumptions, with effective elastic parameters that are assumed to be obtained from 
corresponding experiments. 

Aside from the implication of these effective elastic parameters, we assume that the distributed 
forces, q in vanish. We will focus on the external force that is conjugate to the displacement 
boundary conditions in ([5]) in the rest of this paper. 

Remark 1. The treatment of the kinematics discussed above is in direct analogy with constrained 
formulations in classical, linearized elasticity, such as the plane strain constraint, and the incom- 
pressibility constraint. Consider plane strain: The constraint is assumed to be exactly imposed, 
and we only consider strains in a three-dimensional subspace of the full, six-dimensional space. De- 
note the strains and stresses by and o~ij, respectively. Consider a plane strain problem in which 
the body is loaded by controlling en while letting 022, 012 = and maintaining the plane strain 
constraint, £33, £13, £23 = 0. Importantly, the effective elastic modulus obtained for the an — en 
response differs from the unconstrained case in which, also, £n is controlled while o~ij = for 
i,j 7^ 1. With Young's modulus E and Poisson ratio v for an isotropic material, this effective 
modulus is E/(l — v 2 ) in plane strain, while in the unconstrained case it is E. In plane strain, 
C33 7^ is the Lagrange multiplier enforcing £33 = 0, and ci3,cT23 = are the Lagrange multipliers 
corresponding to £13, £23 = 0, respectively. It will become apparent in Sections [3HH that the only 
loading condition of interest in this paper is exactly analogous to the uniaxial loading discussed in 
this remark in the plane strain context. The identical arguments, point for point, can be made for 
the incompressibility constraint, e\\ + £22 + £33 = 0. 
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3 Force-displacement response of an elastica deforming as a cir- 
cular arc 




ei 



Figure 2: Uncoiling of an elastica shaped as a circular arc. 

Consider an elastica with reference configuration, T, in the form of an arc of a circle with central 
angle 26q and radius R, as shown in Figure [2j The reference positions of points on T are 

C i? (sin O -sin(0 o - S/R)) ) 
X(S) = \ \, Se[0,2R6 } (6) 

( # (cos O -cos(9 -S/R)) J 

In Sections I3.1H3.3I we explore the effect of further assumptions on this deforming elastica. These 
are assumptions of (i) inextensibility, (ii) planar incompressibility of the bounding medium, and 
(iii) stationarity of the strain energy. These cases are not exhaustive. However, the physical 
interpretations and motivations are transparent. 



3.1 The inextensible elastica deforming as a circular arc 

The inextensibility assumption is motivated by the relative stiffness in the arcwise direction of the 
collagen fibril in comparison with the large compliance due to their highly crimped waveform. 

Let the point S = 29qR be displaced by the vector ge\ while the deformed configuration, 7, 
maintains the form of a circular arc without extension, i.e., A = 1. Then, the tip displacement 
is restricted to < g < 2(9q — sin 9q)R. At a given tip displacement, g, the central angle of the 
deformed elastica is 9(g) = 9oR/r(g). The deformed radius, r(g), then satisfies the implicit relation 



r(g) sin 



OqR 



R sin #0 + 2 



(7) 



and the positions of points on 7 are 



x(S) 



r(g) (sin %) - sin [9(g) 


r I cos 9(g) — cos I 9 



(8) 



s 

r(g) 
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From (l6l) and 



d 2 X 

dS 2 



R 



cos (I) 



sin 



(I) 



dx 
dS 





sin (%)-4) 



tfx 
dS 2 



' ^y sin (%)-4) 



k f^y cos (%) 
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(9) 



from which it follows that kq = 1/R, k = 1/r and A = 1. Note that, for g = 2(9q — s'm9o)R, 7 is a 
straight segment of length 29qR along e±. Following (pQ), the strain energy of the elastica can now 
be written as 



20 o R 



W{n-g)= J ^B(K(g)-K Q ) 2 dS, < g < 2(9 - sm9 )R 



(10) 



Note that in addition to W(«; g) being a functional of the field, k, it is a function of the tip 
displacement, g. The force response of the elastica to the tip displacement, g, is 



dW 
dg 



< g < 2(9 -sm9 )R 



(11) 



Like W, the force, /, is a functional of k and a function of g. In what follows, the functional 
character of / is suppressed since its dependence on g is of primary interest. Using k = ||d 2 a;/ds 2 ||, 
and equations ([U]), ([7D and (fTUj) we have, 



fig) 



B9 R 



1 



1 



r{g) 2 \r{g) RJ 9cos9 — sin 



1 



< g < 2(6»o - sin6> ).R. 



(12) 



Equation (fT2|) indicates that f{g) diverges as l/r(g) — > 0. Thus the force in a fully uncoiled, 
inextensible elastica diverges. An extensible elastica, however, develops finite axial tension due to 
stretching along the tangent as it is uncoiled, and will be considered in the next two sections. 

Remark 2. We reiterate that the above approach does not involve a formal solution of (HJ). 
Instead, it is assumed, a priori, that this governing equation is satisfied with the inextensible 
elastica maintaining the circular arc form. 



3.2 The extensible elastica deforming as a circular arc and subject to macro- 
scopic, planar incompressibility 

Collagen fibrils in soft tissues are surrounded by proteoglycan molecules that bind water. At the 
levels of stress that the tissue is subject to, the proteoglycan matrix is nearly incompressible. 
Motivated thus, we consider the waveform of the deforming fibril to be subject to planar incom- 
pressibility as a model for full three-dimensional incompressibility (see Remark 5). The elastica 
deforming as a circular arc in the plane spanned by {ei,e^} is also subject to invariance of the 
area of a circumscribing rectangle, even as the rectangle's aspect ratio varies with tip displacement, 
g. See Figure The conservation of the area, then, leads us to a closed-form expression for the 
height, a, of the current rectangle. 

AgsA ^ 2gy (1 - co.fr) ^ 

2R sin 9q + g 
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2R 2 sin O (1 -cos6» ) 

Figure 3: A circular arc elastica surrounded by a two-dimensionally incompressible medium. 



With this explicit form of a(g), the radius of the deforming elastica can be determined from the 
geometry of Figure [3] as r 2 = (r — a) 2 + (Rsin6o + 5'/2) 2 - This yields 



1 (. , ^, (i? 2 sin0o(l-cos0 o )) 2, 

r (#) = TT-T-T Kf)) + 



(14) 



2«(9) V ~" ' 

The current curvature and stretch A(g) of the circular arc elastica then immediately follow as 

r{g)0{g) 



i 



and A (g) 



R0n 



(15) 



where the implicit relation ((7J) is now replaced by 



0(g) = sin" 



_! i? sin 6»o + g/2 



(16) 



See Figure [3l Inextensibility does not hold: X(g) ^ 1. Using ([IB]) for and A(g), and noting 
that is a function of r(g) and g, we parametrize the strain energy as 



2e R 2 2e fl 





r(9)0(9)X 



1 d5, 



(17) 



where, as previously, PF is a functional of r and a function of g. The tip force, f(g) = dW /dg, is 



f(g) = BR ( - ff 



1 1 

r 



1 \ /3i? 2 sin6» (l - cos6» ) 
V 2a 2 



2i? 2 sin 6 (1 - cos6» ) 



+ 



^-w,9 s / „ ,„ , ( 3R 2 sin 0n(l — cos Oq) a 2 

EAX(X 2 - 1 2sec0 + - tanfl °-\ °-l - ^ 

\ V a rt sm pq ( 1 — cos up 



(18) 

where the derivative formulas dr /dg = 3R sin 6» ( 1 - cos O ) /4a 2 - a 2 /4i? 2 sin 6> ( 1 - cos 6 ) , dO/dg 
sec9/2r — t&n6dr/rdg and dX/dg = (i sec0 + (0 — tan0)dr/dg)/(0o-R) have been incorporated. 
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Remark 3. Note that the stretch, A, defined in (|15p ? is averaged over the reference contour length; 
i.e., over the circular arc with central angle 29q. The use of A in the constitutive relation (|17p implies 
that the axial stiffness, EA, is also homogenized over the reference contour length. 

Remark 4. For this case of the circular arc-shaped elastica in a two-dimensionally incompressible 
medium, it is found that A < 1 for a certain range of macroscopic stretch, A = 1 + g/2Rsm0Q 
(see Section f 5 . 1 j) . The incompressibility constraint causes compression of the elastica. If Q were 
solved with this macroscopic incompressibility constraint it would manifest itself as buckling. We 
have not attempted to follow the buckled shape of the elastica. Instead we have solved for the 
value of #o = ^o cr such that for all 9$ < 9o cr we have dA/dA > 0, ensuring that macroscopic stretch 
translates to microscopic stretch. This root is 9q ci ~ 1-342. Our numerical studies of the circular 
arc-shaped elastica in an incompressible medium (Section 15. ip are restricted to Oq > #o cr . 

Remark 5H The more physically-realistic condition of volumetric incompressibility leads to the 
relation 



from which the results follow in the same manner as outlined in this subsection. The results, how- 
ever, are not qualitatively different between the planar and volumetric incompressibility conditions. 

3.3 The extensible elastica deforming as a circular arc, and relaxing to a sta- 
tionary strain energy configuration 

Referring to Figure [21 and persisting with the stretch, X(g), averaged over the reference contour 
length as in (|15|) 9. the central angle, 9(g), of the deformed elastica as in (|16j) and n(g) = l/r(g), 
the strain energy is still given by (|17|) . 

For a given tip displacement, g, corresponding to an applied force, /, the elastica deforming 
as a circular arc is now assumed to relax to a deformed radius, r(g), at which the strain energy is 
stationary. The motivation comes from the idea that like long chain bio-molecules, a collagen fibril 
also attains an equilibrium state with respect to its configuration, in addition to deforming under 
a tip displacement, g. 

Remark 6. As explained in Section [2.11 the force field q = 0. Therefore, the corresponding work 
term does not enter the stationarity calculations. Furthermore, it is assumed that deformation 
within a subspace of M 3 for displacements, u, manifests in the effective elastic parameters, B and 
E for the elastica deforming as a circular arc. We are interested in stationarity of strain energy 
within this subspace of deformation that the elastica is assumed to explore. That is, the elastica's 
deformation is allowed to vary only over those configurations allowed while maintaining the circular 
arc form for this stationarity calculation. Therefore, the work done by Lagrange multipliers that 
maintain the circular arc shape does not enter the stationarity calculations. With reference to 
Remark 1, an analogous plane strain calculation would seek stationarity of the strain energy within 
the strain subspace {en, £22, £12} while maintaining the plane strain constraint £13, £23, £33 = 0. 
Therefore, the Lagrange multipliers, a 13 , 0*23 , 033 , that enforce these constraints do no work, and 
do not need to be considered in such stationarity calculations. Again, the identical arguments, 

3 It was pointed out to us by a reviewer that the condition in this subsection is properly called a "planar incom- 
pressibility condition" , which we have now done. 




(19) 
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point for point, can be made for a calculation seeking stationarity of the strain energy under the 
incompressibility constraint; i.e., in the subspace defined by E\\ + £22 + £33 = 0. 

The radius, r(g), is the solution of the following stationarity condition: 

<>]] (l UJ + 4EMA(A 2 -l)(0-tan0) = O. (20) 



(21) 



29 BR r> 
or \K r J r 

With the deformed radius thus known, the force is 

dW 



f 



dg 



2EAX (A 2 - 1) sec 9 



since the contribution to / from the derivatives (dW / dr)(dr / dg) vanishes by (|20p . 

Note that the extension of the above results (1121) . (|18j) and (I2ip for half- wavelengths to an 
elastica whose waveform consists of n such half-wavelengths is straightforward. By symmetry, the 
force, /, that results in a total tip displacement of g corresponds to an extension of g/n of each 
half- wavelength, and is obtained from (|12p . ()18f) or (|2ip by substituting g with g/n. 



4 The force-displacement response of a sinusoidal elastica 

For the sinusoidal waveform, the reference configuration of the elastica is defined by two shape 
parameters, the amplitude ao and the half-wave length 1$. See Figure HI As in Section [3] the 




Figure 4: Uncoiling of a sinusoidal elastica. 

elastica lies in the plane spanned by {ei, 63}. The reference positions of points on the elastica can 
be expressed by a single parameter, t, as 



X(t) 




t 


ao sin — t 
V'o 



(22) 



where t E [0, lo] . Analogously, the shape of the deformed elastica can also be formulated in terms 
of a spatial parameter, t, 

< t 





x(t) 




'Ti- 
ll, sin ( —t 



(23) 
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where t € [0, 1] . The deformed half- wavelength, I, is determined by the tip displacement g through 
the relation I := Iq + g. A linear relation therefore exists between t and t 



dt 

t:=tl/l =t(l+g/l ) and — = l/l = (1 + g/l ) ■ 



It then follows that the arguments of the Sine functions in 

-f . 



and 



a {t) := —t 



I 



(24) 

have the same value, i.e. 

(25) 



This geometrical description of the problem allows us to introduce the two main kinematic variables, 
namely the curvature, k, and the stretch, A, as functions of derivatives of X3 and X3 with respect to 
the Lagrangian t and the Eulerian t parameters, respectively. For the planar reference and current 
curves parameterized by t and i, the general curvature formulas given in Section [2] can be simplified 
to the forms 



n (t) :-- 



(l + *3 2 ) 3/2 ' 

The superscript (•)' in ([26]) denotes the derivatives 



*3 



8X3 
dt 



7T 



aoy 1 cos(a(t)) , X% := = -a ^- sin(a(i)) , 



[l + x' 3 y/ 2 



7T 



(26) 







(27) 



7T 



// ._ *9 2 3?3 



7T' 



-a-jj sm(a(t)) . 



b 3 ■- dP 

The local stretch is obtained as A = ds/dS where the infinitesimal arc length measures dS, and 



ds are dS := a/ dX'( + dX| = J 1 + X3 2 dt and ds := y / d~c ! f+~dr| = y 1 + x' 3 2 dt , respectively. 
Then, combining these results with (|24p 9 yields the stretch expression 



A := ds/dS 



1 + 4 z 



l + X 



/ 2 Zq 



(28) 



The basic kinematic variables Ko, n(t, g,a(g)) and X(t, g,a(g)) are thus defined. Note that k and 
A vary pointwise with t and are parametrized by g and a(g). The total energy of the sinusoidal 
elastica in the reference configuration is 

lo lo 
W(k, A; g, a(g)) = J ~B(n(t, g, a(g)) - KQ (t)) 2 J(t)dt + J ^EA(X 2 (t, g, a(g)) - l) 2 J{t)dt (29) 



where J(t) := dS/dt = y 1 + X^ 2 . The tip force f(g), being energy-conjugate to the tip displace- 
ment g, is then given by 



or 



f(g) 



f(g) 



dW(g,a(g)) dW(g,a(g)) 



B(n(t,g,a(g)) - n (t)) 



+ 



+ J EA(X 2 (t,g,a(g))-l) 



dW(g,a(g)) da(g) 
da dg 

J(t)dt 
J(t)dt . 



dK(t,g,a(g)) 



dg 



dX 2 (t,g,a(g)) 



dg 



(30) 



(31) 
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where, as previously, the fact that / and W are functionals of k and A has been suppressed. 

In order to complete the geometric and constitutive description of the deforming sinusoidal 
elastica, a(g) must be obtained for each g. Additional kinematic assumptions are made, as for the 
elastica deforming as a circular arc. In the case of the sinusoidal elastica, the additional kinematic 
assumption determines the current amplitude a. As in Section [3] we consider (i) inextensibility, (ii) 
a planar incompressible bounding medium, and (hi) stationary strain energy. 



4.1 Force-displacement response of an inextensible sinusoidal elastica 

The local inextensibility condition requires that A := ds/dS = 1 and thus 



From ([27]) and {32]) we have, 



A 2 -l = l±4^-l = 0. (32) 



2 2 I 2 -I 2 

tt z cos z (a) 



Clearly, the right hand-side of (|33p can be negative, and is unbounded from below in the limit 
a — > 7r/2. Even for positive values of the right hand-side, i.e. a 2 > 0, a as given by ([33]) varies 
along the elastica. This indicates that the requirement that A = 1 pointwise along the elastica is 
inconsistent with maintenance of the sinusoidal shape. It is worth noting that, even physically, 
it is not clear whether inextensibility should be applied pointwise to model fibrils that are stiff 
to arcwise extension. It may well prove to be more appropriate to consider inextensibility over a 
larger length scale. For this reason, we relax this assumption to a weaker one of conservation of 
the length of the sinusoidal elastica for a given tip displacement g. We continue to require that 
a = a(g) (a function of g only) and therefore is a parameter for the half wavelength. The weak 
inextensibility condition gives 

lo 

j ds = j dS <w J(X(t,g,a)-l)J(t)dt = 0. (34) 

s S 

This condition defines a non-linear residual that can be considered a function of a = a{g) for given 
9- 

lo 

K(a) := j (\(t, g, a) - 1) J(t) dt = (35) 



o 

In order to solve (|35p a standard Newton-Raphson iterative scheme must be employed. Recall that 
this involves the linearization of the residual 1Z{a) about a = a, Lin lZ(a)\ E := 7Z(a) + (dTZ/da)\a(a — 
a) = , and the solution of this equation for a to get a = a — lZ(a) / (d1Z/da)\ a ■ For each given 
value of tip displacement g, this iterative update scheme is repeated until iterates for a converge 
to within a tolerance. Once the value of a is computed, we proceed with the computation of the 
tip force /. To this end, we need the sensitivity of the amplitude a(g) to the tip displacement 
g. It can be calculated by exploiting the implicit form of the residual, now written as 7Z(g,a(g)) 
for a general displacement controlled loading process by writing dlZ(g;a(g))/dg = (dTZ/dg)\ a + 
{d1Z/da)(da/dg) = yielding da/dg = —(d1Z/dg)/{d1Z/da). With this sensitivity in hand, the 
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integrands in (J3TJ) can be computed in a straightforward manner: 



dn 
dg 



Ok 
dg 

dg 



+ 



Ok 
Oa 



9 

d>^ 

da 



do 
dg 



da 




2 ) ^ (1-24 2 ) da 



a 



dg 



(36) 



ad^/ 



4.2 Force-displacement response of a sinusoidal elastica subject to macroscopic, 
planar incompressibility 

As in Section 13.21 we use planar incompressibility as a model for full, three-dimensional incompress- 
ibility. The two-dimensional incompressibility assumption on the surrounding medium (Figure [5]) 
leads to an explicit result for a(g). 




k - 
I 



Figure 5: A sinusoidal elastica surrounded by an incompressible medium. 



a lo = al a{g) = — — = — — . (37) 

I l + g 

Once a(g) is known explicitly in terms of tip displacement g, the derivatives appearing in (|3ip can 
be readily calculated 



dn _ 34' (4 - 1) 3A 2 _ 2/ (1 



and -5T = 12 ,\ . '2 • (38) 



dg I (l + 4 2 )5/2 d 5 / 2 (1 + ^ 2 ) 



4.3 Force-displacement response of a sinusoidal elastica with stationary strain 
energy 

For the sinusoidal elastica, stationarity of strain energy is imposed with respect to the current 
amplitude a(g), i.e. dW(g,a)/da = . This condition defines a non-linear residual 7Z(a) 



11(a) := J (b{k - K )^ + EA(X 2 - J(t)dt , 



(39) 
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which must vanish for a given tip displacement g. As in Section 14.11 Equation f|39|) is solved for a 
by a Newton-Raphson iterative scheme. 

With a being known the tip force, /, can be computed. Since the stationary strain energy 
condition requires that the partial derivative dW(g, a) /da vanishes, the only terms contributing to 
the force will be partial derivatives of kinematic variables with respect to the tip displacement g. 
That is, it suffices to compute the integrand terms 

8k _ gg (4 2 - 2) d\ 2 _ 2 1 1 

^"T(i + x ^ )5/ 2 and -dj-q (TT^y (40) 

Remark 7: Consider the limiting case in which the tip displacement g is much larger than the 
reference half-wavelength Iq, i.e. g/lo S> 1, and I/Iq & g/l Q 1 . The elastica tends toward the 
limiting shape of a straight segment along e\. Owing to the flat shape of the elastica, its spatial 
slope x' 3 and the curvature xg are small. This implies that the contribution from the bending term 
to the tip force in both (|38P i and ()4Q[) i will be negligible in comparison with the contribution 
from the axial extension. Furthermore, the vanishing term x' 3 in ([38])2 for large values of the tip 
displacement g makes the force terms in (|38p ? and (|40p 9 tend toward each other. Then, provided 
the bending stiffness is not much larger than the axial stiffness, the tip force values for the planar 
incompressible bounding medium and stationary energy cases approach a common value at large 
tip displacement values. This is reflected in Figures [12] and [13j 

Remark 8: Computations with all the models of the sinusoidal elastica discussed in the preceding 
Sections 14. If - 14.31 require an efficient numerical integration tool both for computing the tip force, 
/, and for carrying out the Newton-Raphson iterations. For this purpose, we have employed a 
set of F77 subroutines, the so-called dcuhre, providing a double precision integration tool based 
on adaptive division of the integration domain into subre gions. For deta i ls of t he theory and the 
implementation of the algorithm, the reader is referred to Berntsen et al. ( 199 ll ). 



5 Comparison of shape, kinematic and stationary energy assump- 
tions; validation 

We now turn to a comparative study of the force-displacement response of the circular arc and 
sinusoidal elasticas, with the additional goal of gaining insight to matches between the models and 
experimental data. 

5.1 The force-displacement response of circular-arc and sinusoidal elasticas sub- 
jected to different kinematic assumptions 

We first consider the force-stretch behavior of the circular-arc elastica subject to the two additional 
kinematic assumptions and the stationary strain energy assumption. To this end, we first relate 
the micro-tip displacement, g, and the macro-stretch, A. The displacement between the ends of 
the elastica is assumed to be dictated by macroscopic deformation in an affine manner. That is, the 
macro-stretch A is related to the tip displacement g by A := 1 + g/{2R sin 6q) (see Figure [2]). □ In 
the studies to follow, the macro stretch will be used as the primary deformation variable controlling 
the force. 

4 For tissues with transverse isotropy, where the collagen fibrils (elasticas) are characterized by end-to-end vectors 
that are highly aligned, affinity of deformation is a good assumption. The alternative, fibril slippage, will be treated 
in a separate paper. 
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a) 0o b) \ 

Figure 6: Circular-arc, inextensible elastica. a) Variation of the maximum macro-stretch Ah ee i = 
9q/ sin#o with the initial angle 9q. b) Shift of the maximum value of the attainable macro-stretch 
Aheei for selected values of the initial angle 6q = | , ^ , ^ , and for B = 1, R = 1. 

For the circular-arc, inextensible elastica, according to (|12|). the tip force / diverges as the tip 
displacement g approaches the value <7 max = 2-R(0o — sin#o)- This implies that the maximum value 
of macro-stretch is Aheei = 00 / sin0o , as shown also in Figures [6^ and [6b. The geometric parameter 
00 thus has an unambiguous physical effect. The ordinates of Figure can be obtained as the 
locking stretch for each value of 9q as shown more transparently in Figure [Ub for the /—A response 
parametrized by 8q. For the circular-arc, inextensible elastica, the shape of the f—X curve depends 
only on the initial radius R and the bending stiffness B. Figures [T^t and [7b demonstrate that the 
larger the initial radius R or the smaller the bending stiffness B, the sharper is the transition to 
divergence in /. Although the sharpness of the transition can be tuned, the value of the locking 
stretch remains Aheei = &o/ sin 9q , and the response beyond the heel region is asymptotic to a 
vertical line at Aheei- Of course, this divergent f—X response is non-physical, and considerably 
limits the ability to match experiments on different collagenous materials possessing distinct, and 
non-divergent, responses in the post-heel region. 

The variation of the micro-stretch, A, with macro-stretch, A, for a circular-arc elastica embedded 
in a planar incompressible medium is depicted in Figure [8^. We draw attention to the compression of 
the elastica for a regime of deformation characterized by small values of A, and discussed in Remark 
4. Following the approach outlined there we have solved for 9q = 0o cr such that for all 0q < 0o cr 
we have dA/dA > 0. Explicitly we_have A = 0(X)r(X)/(R9 o ) with r(A) = R(l - cos(0 o ))/(2A) + 
A 3 J Rsin 2 (6» )/(2(l -cos(0 o )) and 0(A) = sin^ 1 (Ai?sin0 o /r(A)) . Setting the derivative dA/dA|x = i = 
and solving the resulting equation for 8q, we obtained the maximum: 0q < 0q ci ~ 1.342. As 
clearly shown in Figure [8^,, for the values of the initial angle 9q < 0q ci the compressive regime of 
micro-stretch is avoided. For this reason our subsequent investigations are restricted to 8q < 0o cr 
for the circular-arc elastica embedded in a planar incompressible medium. 

The sensitivity of the tip force, /, to the axial stiffness, EA, for B = 1, R = 1 and 8q = tt/3 is 
depicted in Figure [8b for the circular-arc elastica surrounded by a planar incompressible medium. 
A decrease in the axial stiffness translates, as expected, to a decrease of the slope. Clearly, the 
axial stiffness dominates the slope of the f—X curve. In contrast to the inextensible case we observe 
neither a long toe region nor a distinguishable heel region. Larger values of A are attainable. As 
the axial stiffness is increased, however, the f—X curve enters the high (but still increasing) slope 
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Figure 7: Circular-arc. inextensible elastica. Sensitivity analysis of the f — X curve to a) the 
variation of initial radius R G [0.4, 1] and b) the variation of bending stiffness B € [1, 10] for 
6 = ir/2. 



regime about Aheel = Go/ sin#o without exhibiting much of a toe region. This, of course, limits the 
use of this model. 

The material parameters used for the circular-arc elastica with stationarity of strain energy 
were the same as for the planar incompressible medium case. As in the planar incompressible case 
a decrease in axial stiffness has a strong, depressing influence on slope of the f—X response. The 
elastica can be extended to A > Oo/smOo • With an increase in axial stiffness, the f—X curve 
approaches the behavior of the inextensible circular-arc elastica (see Figure [9]). The circular-arc 
elastica attaining a stationary strain energy possesses a number of favorable properties: The toe 
region exists and its slope can be tuned by the bending stiffness. The location of the heel region is 
uniquely determined by the initial angle 6q as Ah ee i = Go/ sin#o- The slope of the post-heel region 
can be adjusted by the axial stiffness EA as shown in Figure [9J The stationary strain energy 
assumption with clearly identifiable parameters thus serves as a promising model to match with 
experimental data. 

In Figure [10] we compare all three cases of the circular-arc elastica. Distinct values EA = 
34, 67, 100 are assigned to the axial modulus of the planar incompressible and stationary energy 
elasticas. The f—X curve of the stationary energy case approaches the inextensible one by "rotating" 
about the heel just below Aheel- However, this occurs with no discernible difference in the curves 
for A values smaller than the heel. In case of the elastica surrounded by a planar incompressible 
medium, however, the stiffening in the f—X behavior is different. Owing to larger values of micro- 
stretch in the initial stages, the location of the heel shifts to smaller values of A. 

In the foregoing parameter study, we solely considered circular-arc elasticas with two kinematic 
assumptions and the stationary strain energy assumption. In what follows, we present an analogous 
parameter sensitivity study for the sinusoidal geometry. In contrast to the circular-arc elastica, the 
reference shape of a sinusoidal elastica is governed by two parameters: the amplitude ao and 
the half-wave length Iq (see Figure fH). The ratio an I In, however, cannot be arbitrarily chosen. 



According to the results reported bv lDale et al.l ()1972j ). this ratio is limited to values smaller than 



0.1. Accounting for this fact in the studies to follow the ratio has been chosen as clq/Iq < 0.2 , which 
will allow us to consider values slightly larger than the experimental observations. The macro- 
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Figure 8: Circular-arc elastica surrounded by a planar incompressible medium, a) Dependence of 
the micro-stretch A on the macro-stretch A and the initial angle 9q. b) Sensitivity analysis of the 
/—A curve to the variation of the axial stiffness EA £ [1, 100] . curve. 



stretch, A, remains the primary deformation measure, and is now related to the tip displacement, 
9, by A = l + g/l . 

First, we consider a sinusoidal elastica with the additional global inextensibility assumption 
given in ([34]) . In Figure [TTa the influence of the ratio ao/lo G [0.05,0.2] on the /—A curve is 
depicted while keeping the material parameters fixed at B = 1 and EA = 1. This ratio proves 
crucial in determining the value of stretch at which the heel occurs. The higher the ratio ao/lo, the 
longer the toe region preceding the heel. In other words, this parameter determines the value of 
A where the influence of the bending mechanism starts to diminish and the axial extension begins 
to govern the /—A curve. In order to demonstrate the sensitivity of the /—A curve to the bending 
stiffness, the ratio of bending stiffness to axial stiffness, B/EA, is varied from 1 to 4 (Figure [Tib). 
An increase in the ratio B/EA scales the curve's ordinates (/-values), and therefore the transition 
in the heel region becomes more gradual. However, the value of the locking stretch is not influenced 
by the changes in the ratio B/EA. 

In the last two cases we consider the planar incompressible and stationary energy sinusoidal 
elasticas. Figures [T2a and 113a present the influence of the change in ratio ao/lo on the /—A curves 
of the respective cases. Like the inextensible case the ratio ao/lo is varied within the interval 
[0.05,0.2] while the value of the ratio EA/B is kept fixed at 30. Clearly, the /—A curves for the 
planar incompressible and stationary energy cases do not exhibit a sharp transition to stiffening 
behavior. This is in contrast with the inextensible case in Figure [TIT Variation of the ratio ao/lo 
does not cause significant change in the shape of the curves. 

The sensitivity of the /—A curves for the separate cases to changes in material parameters 
EA and B is presented in Figures 112b and 113b . respectively. The ratio EA/B varies in the 
range [30,300]. The axial stiffening is clearly reflected in the curves. No striking shape change is 
observed. We draw attention to the fact that the /—A curves in Figures [T2l and [T3l for the planar 
incompressible and stationary energy cases of the sinusoidal elastica are quite similar. The reasons 
for this similarity have been already outlined in Remark 8. 
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Figure 9: Circular-arc elastica with stationary energy. Sensitivity analysis of the /—A curve to the 
variation of the axial stiffness EA S [1, 100] . 



5.2 Comparison with experiment 

In the preceding section the sensitivities of the / — A curves to geometric and material param- 
eters have been discussed for both the circular-arc and the sinusoidal elasticas subjected to two 
additional kinematic assumptions, and the sta tionary strain energy assum ption. In this section we 
carry out a comparison with data reported by Freed and Doehring (2005) (Figure [H]). These data 



correspond to uniaxial extension experiments on five chordae tendineae from procine mitral valves. 
They demonstrate a long toe region relative to the maximum stretch in each experiment. At the 
heel (Aheel ~ 1-13), the nominal stress-stretch curve stiffens sharply to a larger slope. From the 
results in Figures I6HT31 of the preceding parameter study, it is apparent that this behavior can only 
be captured either by the inextensible sinusoidal elastica, or the circular-arc elastica attaining a 
stationary energy state. Figure [TH compares the experiment with these two models with the mate- 
rial parameters given in the caption. Both the sinusoidal and the circular-arc models successfully 
match the data in the toe region. The inextensible sinusoidal model can also predict the upturning 
region, but its stiffness beyond the heel region rapidly diverges and fails to match the experimental 
results. In the case of the circular-arc model, the value of 9q can be analytically determined from 
the macro-stretch value at the heel. We solve for 0q such that 6q/ sin#o = Ah cc i = 1-13. This gives 
the initial angle #o ~ 4-7r/15, and the ratio of the axial and bending stiffness EA/ B can be tuned 
to match the slopes of the regions just preceding and succeeding the heel region. The initial radius 
R is varied to match the sharpness of the slope change at the heel region. The comparison of the 
stationary strain energy circular-arc elastica and the experimental data clearly illustrates that the 
proposed model quantitatively captures the experimental data with just a few parameters: 9o, R, B 
and EA, all of which are very well-motivated physically. Clearly, other such experimental data can 
be matched without difficulty. 



19 




Figure 10: Comparison of circular-arc elasticas subjected to different constraints. In the planar 
incompressible and stationary strain energy cases, three different values are assigned to the axial 
stiffness EA = 34, 67, 100, corresponding to increasingly stiff /—A response. 



6 Macroscopic material model incorporating the elastica 

6.1 Continuum strain energy density function at the macroscale 

The contribution to the overall strain energy density function due to the collagen fibrils embedded 
in a nearly incompressible viscous medium is obtained by summing up the free energies of individual 
elastica- like fibrils, 

*coi = -^rW(g) . (41) 

With the unit vector e denoting the average orientation of collagen fibrils, the macroscopic stretch 
in this direction is obtained by A = \Fe\, where F is the deformation gradient tensor. In the 
context of anisotropic elasticity, especially transverse isotropy, it is common to define structural 
tensors M := e <g> e for the construction of strain energy density functions formulated in terms of 
additional invariants. The derivatives of these invariants are then used as tensor generators in the 
stress response functions. In the present case, I4 := C : M = A 2 is the relevant invariant (C being 
the right Cauchy-Green tensor). We continue to use an affine relation between the macro-stretch 
and tip displacement of fibrils, i.e. A = 1 + g/lo , where 1$ is the half wavelength of a sinusoidal 
elastica, and Iq = 2i?sin#o for a circular-arc elastica. The cross-sectional area of the tissue that 
contains N such fibrils is Aq. With this relation in hand, the contribution to the total second 
Piola-Kirchhoff stress tensor due to the stretching of collagen fibrils, S co1 = 2d^ co \/dC, can be 
obtained as 

5 coi = N /(fl)M M ^ 
Aq 

where the results dg/dX = Iq , 2d\/dl4 = 1/A , dl^/dC = M and the definition f(g) := dW/dg 
have been used. Then, the nominal stress tensor P co1 readily follows from P co1 = FS co1 : 

pcol = N [(3le®e (43) 
Aq 

where Fe = Ae and lei = 1 . 
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Figure 11: Sinusoidal inextensible elastica. Comparison of the f—X curves for globally inextensible 
sinusoidal elasticas having different a) ao/^o and b) B/EA ratios. 



500 



400 



300 



■aoAo = 0.2 
aoAo = 0.14 
aoAo = 0.08 

■ ao/Jo = 0.02 

SA/B = 30 



200 



100 




500 



400 



300 



200 



100 



( 


i / 


' : / 
/ ■ / 
i ■ ' 


/ — EA/B = 30 




/ -- EA/B = 120 




/ ■■■ EA/B = 210 




/ — EA/B = 300 




ao/io =0.1 



a) 



6) 



1.5 



2.5 



Figure 12: Sinusoidal elastica surrounded by a planar incompressible medium. Comparison of the 
f—X curves for different a) clq/Iq and b) EA/B ratios. 



We now turn our attention to the convexity of the strain energy density = ^coi(^4)) i n 
order to have a basic understanding of its stability properties. The convexity condition demands 
the positive definiteness of the first elasticity tensor A co1 

H : A co1 : H > VJET € M 3x3 and A co1 := , (44) 

where M 3x3 is the space of second-order tensors in M 3 . 

The explicit form of the first elasticity tensor A co1 can be obtained by the chain rule as 

A ■■=^W0F + ^dF^dF (45) 
where the superscript (■)' denotes the derivatives with respect to I4. The quadratic product of A co1 



21 




Figure 13: Sinusoidal elastica deforming by attaining a stationary strain energy state. Comparison 
of the f—X curves for different a) qq/Iq and b) EA/B ratios. 



with H involves the terms 

H; ^H? :H = 2| ' He| ' 2a0 - (§:") 2 = (2He.Fe)>>0. («) 

Note that both terms are non-negative. The local convexity condition (|44|) of the free energy ^> co \ 
reduces to 

2*' col \\He\\ 2 + 4^' ol (He ■ Fef > . (47) 

Based on (|4U|) . non-negativity of both ^' col and is sufficient to fulfill the convexity condition 
(|4T1) . though not necessary. The explicit forms of the derivatives are = Nf(g)/(2Ao\) and 
^" ol = Nlo(f'(g)lo + gf'(g) — f (g)) / (4AqX 2 (g + Iq)). If we assume that collagen fibrils can carry 
only tensile loads, the positiveness of W j is satisfied identically for f(g) > . Furthermore, the 
convexity of ^ co \ with respect to g ensures that f'(g) > . Thus, it is now sufficient to show that 
the term gf'(g) — f(g) > in V'L . This condition can be obtained starting from the convexity 
condition for f(g) with respect to g, i.e. f"(g) > . For positive values of g, we have gf"(g) > . 
Integration of gf'(g) > by parts yields f° gf'(g) = gf'(g)\ 9 - J S /'(<?) >J) • For /(0) = 0, we 
obtain the sought form gf'(g) — f(g) > . Therefore, convexity of both the W and f(g) = dW /dg 
guarantees the local convexity of the macroscopic free energy function \I/ co i . 

When O co i is combined by rule-of-mixtures with the strain energy density function of the sur- 
rounding matrix medium, the above results completely characterize the influence upon the convexity 
of the overall composite, leaving open only the question of convexity of the matrix material. How- 
ever, as pointed out by a reviewer, the actual interaction between the collagen fibrils and matrix 
involves shearing of the matrix and consideration of the rate of decay of shear fields with distance 
from a loaded fibril. There is then the possibility of more complex interaction between strain energy 
of the matrix and of the elastica- like fibrils. The convexity arguments presented here will then have 
to be refined by further considerations of fibril-matrix interactions. 
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Figure 14: Simulations of experimental data by the inextensible sinusoidal elastica (aoAo = 
0.245 ,B/EA = 25 mm -1 ) and the circular-arc elastica attaining a stationary strain energy state 
(R = 0.013 mm , 0q = A-k/15 , EA/B = 7 x 10 6 mm -2 ). In the elastica models R, clq and 1$ were in 
mm B in N/mm and EA in MPa/mm. 

7 Closing remarks 

The primary aim of this paper is a discussion of the characteristic soft tissue response in the context 
of the elastica-like mechanical behavior of slender fibrillar structures in these tissues. The models 
are applicable to tendons, skin and the passive response of muscle. While entropic elasticity-based 
models can also model this characteristic soft tissue response — especially the locking behavior — 
there are strong physical and physiological reasons to surmise that this is the wrong approach 
to adopt. A direct solution of the shape and force in a deforming elastica requires the solution 
of a ("highly") nonlinear, fourth-order partial differential equation. The simplification used here 
is that judiciously-chosen additional assumptions on the kinematics and on the energy state can 
lead to force-deformation response functions for the elastica. This is the central thesis advanced 
in this paper. Beyond this, the paper is concerned with an enumeration of two families of shapes 
(circular arcs and sinusoidal half-periods) of the deforming elastica, and three possible additional 
assumptions: inextensibility, macroscopic planar incompressibility, and stationarity of strain energy. 
The motivations for each of these additional assumptions are well-founded in a physical sense. 
Their suitability in matching a set of experimental force-deformation curves has been examined. 
On the basis of the current limitation to elastic effects, it emerges that the elastica deforming as 
a circular arc, and maintaining itself in a state of stationary strain energy in each configuration 
(parametrized by overall elongation) can resolve the experimental data to a high degree of precision!! 
The parameters used are the two stiffnesses — bending and axial, and two geometric parameters that 
determine the shape of the undeformed elastica. These can be easily determined from mechanical 
experiments and micrographs, and compared with the values obtained for the best fit. Such an 
exercise would be a strong validation of these models. We note that the circular-arc elastica 
with stationary energy matches the experimental data very well in Figure [14] with initial radius 
R = 0.013 mm (13 /im), and initial central angle 26q = 87r/15 (96°). This gives a wavelength of 
4i?sin#o = 38.64 /im which seems very reasonable, given that collagen fibrils are typically found to 

5 We have not demonstrated quantitative error measures since no significant physical insight is gained by doing so. 
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have w avelengths between 10 and 50 //m as in IScreen et al. I (l2004h and iProvenzano and Vanderbvl 



(2006). The importance of matrix shear lag and its influence on convexity, inelastic effects such as 



the viscous friction as collagen fibrils move relative to the surrounding proteoglycans, viscoelasticity 
of the collagen fibrils themselves and proteoglycans, and slippage of fibrils under larger forces, must 
not be overlooked, however. 

It should also be quite clear, that the development here has complete relevance for many types of 
slender filamentous structures, from carbon nanotubes, through underwater cables to oil pipelines. 
The class of continuum strain energy density functions so developed in Section [6] is applicable to 
any composite consisting of mainly unidirectional, elastica-like reinforcing fibers in a matrix. 
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